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We consider a simple integrable model of a spin chain exhibiting the Magnetoelectric Effect 
(MEE). Starting from the periodic 5 = 1/2 XXZ-chain with Dzyaloshinskii-Moriya terms, which 
we consider as a local electric polarization in the spirit of the Katsura-Nagaosa-Baladsky (KNB) 
mechanism, we perform the mapping onto the conventional XXZ-chain with twisted boundary condi- 
tions. Using the techniques of Quantum Transfer Matrix (QTM) and Non-Linear Integral Equations 
(NLIE) we obtain the magnetization, electric polarization and magnetoelectric tensor as functions of 
magnetic and electric field for arbitrary temperatures. We investigate these dependencies as well as 
the thermal behavior of the above mentioned physical quantities, especially in the low-temperature 
regime. We found several regimes of polarization. Adjusting the magnetic field one can switch the 
system from one regime to another. The features of the critical properties connected with the MEE 
are also illustrated. 
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I. INTRODUCTION 

The strong interest in multiferroic materials (materials 
exhibiting simultaneously several primary ferroic order 
parameters), which occurred recently, ^""^ also triggered 
intensive investigations in the field of magnetoelectric ef- 
fects (MEE) and magnetoelectric materials.'* In general, 
MEE denominates the mutual influence of the electric 
and magnetic properties in matter. In its most promi- 
nent form, MEE can be defined as the magnetic field 
dependence of (ferro)electric polarization and the elec- 
tric field dependence of magnetization. In linear ap- 
proximation, one can introduce the linear magnetoelec- 
tric tensor (isothermic), which quantifies the induction 
of (ferro)electric polarization by a magnetic field and of 
magnetization by an electric field: 



T, E 



(1.1) 



T, H 



where T is the temperature, i,j = {x,y,z), P (M) 
the (ferro)electric polarization (magnetization) and E 
(H) the electric (magnetic) field, respectively. The ef- 
ficient control of magnetic properties of solids by means 
of the electric field has a lot of potential applications, 
for instance, in spintronics.^ At the moment, the num- 
ber of known materials exhibiting the MEE exceeds one 
hundred.^ One of the microscopic mechanisms of the 
MEE was theoretically described by Katsura, Nagaosa, 
and Baladsky in 2005 (KNB-mechanism)^ and connects 
the appearance of the local (ferro)electric polarization 
with the non-coUinear ordering of the neighboring spins, 
i. e. with the spin current 



where is the unit vector connecting two neighbor- 
ing spins Si and Sj. The general phenomenological ar- 
guments of Ginzburg-Landau theory^ as well as sym- 
metry arguments^" showed that this coupling between 
(ferro)electric polarization and magnetization always ex- 
ists in certain classes of materials independent of the 
crystal symmetry. The MEE has been successfully de- 
scribed within the KNB-mechanism for many classes of 
ferroelectric materials. -'^^^^'^ In turn, the discovery of the 
deep relationship between magnetic structure and ferro- 
electricity in several spin-chain materials stimulated fur- 
ther research of multiferroics and paved a way between 
the physics of multiferroics and quantum spin systems. 
The simplest system exhibiting the MEE by means of 
KNB-mechanism, considered so far, is the S = 1/2 Ji- 
J2-spin chain, which is believed to be a more or less ad- 
equate model for several multiferroic materials such as 
LiCuaOa,*^"^^ LiCuV04,*^"2° CuCls,^! and others. The 
expression for the electric polarization becomes particu- 
larly simple in the case of a one-dimensional linear chain, 
say, in x-direction. In this case, all e-y = ei^i^i = e^,, 
and one gets for the electric polarization according to 
the KNB-mechanism: 



P. = 0, 



n=l 
1 ^ 



(1.3a) 
(1.3b) 

(1.3c) 



Gij X (Sj X Sj) , 



(1.2) 



where 7 is a material-dependent constant. Thus, con- 
sidering one of the simplest lattice spin models exhibit- 
ing some interplay between magnetic and (ferro)electric 
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properties, one arrives at the Hamiltonian of a J1-J2- 
chain in a magnetic field with additional Dzyaloshinskii- 
Moriya (DM) interaction. ^^"^^ Supposing the magnetic 
field to be in z-direction and the electric field in y- 
direction, one gets: 



N 



N 



H — Ji S„S„+i + J2 S„S„+2 (1-4) 

n— 1 n—1 

N N 

n—1 n—1 

Although the frustrated S = 1/2 spin chain with ferro- 
magnetic nearest-neighbor and antifcrromagnctic ncxt- 
to-nearest-ncighbor interactions has recently received 
considerable amount of attention, ^-^"^^ only few theoreti- 
cal (mainly numerical) works have been devoted to MEE- 
connected issues and to the interaction with an electric 
field by means of DM-terms. The full Hamiltonian (1.4) 
is not integrable, whereas our aim is to develop an ex- 
act analysis of the MEE. Therefore, we look at a slightly 
simpler model with only nearest-neighbor interactions. 
Thus, we consider the 5=1/2 XXZ-chain with DM- 
terms. We also suppose the potential ability of the sys- 
tem to generate an electric polarization due to the KNB- 
mcchanism. The integrability of the 5*= 1/2 XXZ-chain 
with DM-terms is based on a gauge transformation to 
the ordinary XXZ-chain. For the latter we apply the 
method of the quantum transfer matrix (QTM) leading 
to non- linear integral equations (NLIE)^^"^° which deter- 
mine the exact free energy of the system, and thus, we 
are able to obtain exact thermodynamic functions. 



II. THE MODEL 

A. Hamiltonian and MEE parameters 

Let us consider the 5 = 1/2 XXZ-chain in a longitudi- 
nal magnetic field and with Dzyaloshinskii-Moriya inter- 
action: 

N 



one obtains the Hamiltonian of the ordinary XXZ-chain 
with certain changes in J and A and with twisted bound- 
ary conditions. In other words, the gauge transformation 
corresponding to the site-dependent rotations of spins is 
equivalent to the following canonical transformation of 
the initial Hamiltonian:^^^^^ 



n = unu 

N 



(2.4) 
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N 



Z^ = exp|-z0^j5| \ , 
where 



n=l 

(2.5) 
(2.6) 



As the Hamiltonian of Eq. (2.4) is integrable and the 
methods of QTM and NLIE provide a very efficient for- 
malism to calculate finite-temperature properties of the 
model, we are able to describe the MEE by the influence 
of magnetic and electric fields on thermodynamic func- 
tions. We are especially interested in the magnetization 
and the electric polarization which are defined as follows: 



n=l 
1 ^ 

^ ~ ~fj E {SnSn+1 ~ S^L^n+l) i (2-7b) 



Here, the angle brackets denote the thermal average with 
respect to the statistical operator p = e~'^^'^ /Z. These 
quantities can be derived from the canonical partition 
function Z = tr(e^^/-'") of the system. Taking the con- 
nection between the Hamiltonians (2.1) and (2.4) into 
account, one can write down the following relations: 



n = jY^ {S^S^i^, + SlSl^, -f A5,15,l+0 (2.1) ^„^(y^ H,E) = -(^ 



T.E 



N 



N 



J^'^^ {.SnSn+1 'S'n'S'n+l) ^ '^n ■ 



n=l 



n=l 



Here, periodic boundary conditions arc assumed, S'^j^^ = 
5f , and the material-dependent constants 7 and g^,B are 
absorbed into the electric field strength JE = i?y7 and 
the magnetic field strength H = gp^sHz^ respectively. 
The operators 5" obey the standard SU (2)-algebra, 



P{T,H,E) j\dEj^H JdE\dJ^ 




, (2.8a) 



E 



G^^(T,i/,J,A)^ 



a{T,H,E) 



_lfdM\ f dP 



(2.2) 



J \ dE J u \ dH / jt; 



T,H,.I-A 



(2.8b) 



(2.8c) 



Implementing the rotation of all spins about the z-axis 
by angles proportional to the number of lattice sites, 



+inch 



tan~^(£') 



(2.3) 



where / and / are the free energies per lattice site for the 
conventional XXZ-chain (Eq. (2.4)) with i?-depcndent 
coupling constant and anisotropy (Eqs. (2.6)) and for the 
XXZ-chain with DM- interaction (Eq. (2.1)), respectively. 
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Note that the derivative in Eq. (2.8b) with respect to J 
has to be performed at constant value of the product J- A. 
Alternatively, we may apply the canonical transformation 
(2.5) to the definition (2.7b) to obtain the result (2.8b), 
where G^** is the nearest-neighbor transverse correlation 
function: 



N 



N ^ 

ri=l 



(2.9) 



Here, angle brackets denote thermal averaging with re- 
spect to the statistical operator p = e~'^/-^/tr(e~^/-^) . 



B. Phase diagram 

The correspondence between 5* = 1/2 XXZ-chains with 
and without DM-interaction allows one to map the well 
known ground-state phase diagram of the XXZ-chain 
onto the (E^ i?/J)-plane (see Fig. 1). The region between 
the two curves corresponds to the massless polarized (an- 
tiferromagnetic) ground states, and the region below the 
lower curve corresponds to the Neel ordering of the spins. 
The following equations describe the iJ-depcndence of the 
upper and lower critical fields: 



i/„/j = A + v^TT^, 



k— — oo 



ch (fcry) ' 



(2.10a) 
(2.10b) 



III. RESULTS 



Free-fermion case 



As known, in case of zero anisotropy A = 0, the S = 
1/2 XXZ-chain reduces to the much simpler XX-chain 
which, in turn, can be solved exactly by mapping it onto 
a free-fermion system by means of the Jordan- Wigner 
transformation.^^ The XX-chain with DM-term is also 
exactly solvable within this context. Thus, starting 
with the Hamiltonian (2.1) with A = and performing 
a Jordan- Wigner transformation to free fermions: 



(3.1a) 



Sn 


n-1 

= n(i 


- 2c+c,) c+ 










n-1 


- 2c+cj) c„ 




i=i 




qZ 


— Cfi ^11 ~ 


1/2, 



(3.1b) 
(3.1c) 



where and Cn are creation and annihilation oper- 
ators of spinless fermions obeying the standard anti- 
commutation relations, {c„,Cm} = 0, {c^,Cm} = 5nm, 
one arrives at the following free-fermion Hamiltonian: 



NH 



N 



DM 



(3.2) 



N 



where ch rj = A = A/y/l + E^ is an appropriate 
parametrization of the anisotropy A > 1, i. e. i? < = 
VA^ — 1. For E > El the electric field is too large for 
Neel ordering. 
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FIG. 1: (Color online) Ground-state phase diagram of the 
5 = 1/2 XXZ-chain (2.1) for A = 4.5 in the (S, ii"/ J)-plane. 
For i5 = the phase transitions occur at Hi/ J = 2.637. . . 
and Hu/ J = A -I- 1 = 5.5. For E > Ei there does not exist 
any Neel ordering. 



Here, as usual, periodic boundary conditions are as- 
sumed. Then, the Fourier transform diagonalizes the 
Hamiltonian and yields 

^xx-DM = E^*^ ("^fe^*^ ^ -^/^) ' {^."ia) 

k 

ek = -H + J{cosk + Esink) . (3.3b) 

For the free-fermion picture of the model under con- 
sideration, one can easily obtain all thermodynamic func- 
tions in the thermodynamic limit in the form of integrals 
over the one-dimensional Brillouin zone, k G [— tt, tt] (see 
for example Ref. [50]). So, one gets for the free energy / 
per spin, the magnetization M, the polarization P and 
the magnetoelectric tensor a: 



f{T,H,E)^l- 1" (^+Tlog(nfc))d/c, 
M{T,H,E) = ^ f {nk-l/2)dk, 

27r J-yr 

{ilk — 1/2) sin/c d/c , 



p(r, ij, E) 

a(T, H, E) 



1 
1 



27rr 



(3.4a) 
(3.4b) 
(3.4c) 

7ife(l — nfc) sinfcdfc , (3.4d) 
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FIG. 2: (Color online) Magnetization M and polarization P 
as functions of the magnetic field H for anisotropy A = 0, 
temperature T/J = 0.01 and different values of the electric 
field E = 0,1,2,3. 



where Uk = 1/(1 + e'^''^^) are the occupation numbers of 
spinless ferniions. 

Figure 2 shows the dependence of the magnetization M 
and the polarization P on the magnetic field H for dif- 
ferent values of the electric field E in units of J and for 
very low temperature. The phase transitions at the corre- 
sponding values of iJ = Hu = vl + E^ can be observed. 
All critical exponents arc 1/2. The curves in Fig. 2 are 
obtained by the method described in the next section and 
in App. A, see for instance Eqs. (3.5) and (3.7). Simi- 
larly, direct numerical calculation of the same quantities 
by computing the integrals on the r. h. s. of Eqs. (3.4) 
produce data with deviations smaller than 10~^. 



B. General case 

It is convenient to express the magnetization and the 
polarization in terms of static short-range correlation 
functions of the XXZ-chain. Using Eqs. (2.7a), (2.8b), 
and (2.9) as well as the translational invariance of the 
Hamiltonian (2.4) one obtains 



M{T,H,E) = {S! 
p(r, H, E) : 



E 



(3.5a) 
(3.5b) 



Due to recent progress in the theory of intcgrablc mod- 
els static short-range correlation functions of the XXZ- 
chain can be calculated exactly. It can be shown that all 
static correlation functions of the one-dimensional XXZ- 
model can be expressed as polynomials in the derivatives 
of three functions^^ ip, uj, and uj' which are determined by 
solutions of certain numerically well behaved linear and 



non-linear integral equations. In App. A we give the 
definitions of these functions in the massive case A > 1.^^ 
In the critical case < A < 1 the definitions are quite 
similar. 

Employing the short-hand notations 



= 9^<y5(a;)U=o , 
/(m,«) = d^dyf{x,y)\^=y=o for / = a;,w' , 



we obtain^'^ 



J{0 



,0) ch77w'(i 



0) 



i sh?7 



8ry 



(3.6a) 
(3.6b) 



(3.7a) 
(3.7b) 



where the parameter rj is defined by A = ch 77 for A > 1 , 
i. e. E < \J A^ — 1. Similar expressions can be obtained 
for the case < A < 1,^4 i. e. £; > VA^ - 1. 

Combining Eqs. (3.7), (3.5), and (2.8c) the magneti- 
zation, the polarization and thereby the magnetoelectric 
tensor can be calculated exactly for all temperatures over 
the whole range of the phase diagram Fig. 1. 

By solving the linear and non-linear integral equations 
of Appendix A, the MEE parameters M, P and a of 
Eqs. (3.5a), (3.5b) and (2.8c), respectively, can be calcu- 
lated with high numerical accuracy. 

Figure 3 shows these three quantities as functions of 
the magnetic field for very low temperatures and differ- 
ent values of the electric field. The phase transitions at 
the corresponding values of the magnetic field (see Fig. 1) 
can be easily observed. Interestingly, in the Tomonaga- 
Luttingcr liquid phase, the electric polarization P is a 
non-monotonic function of the magnetic field H whereas 
it is constant, P 0, in the Neel-ordered phase and van- 
ishing in the fully polarized ferromagnetic phase. Here, 
the critical exponents characterizing the magnetic field 
dependence of the electric polarization in the vicinity 
of critical values of the magnetic field, H = Hi{E) and 
H = Hu{E), are all equal to 1/2, 



P^{H^ Hi) 
P ^ (H^ - H] 



1/2 
,1/2 



(3.8a) 
(3.8b) 



This can be also observed by the van Hove-like singu- 
larities of the magnetoelectric tensor a{H). We evalu- 
ated the magnetoelectric tensor a(r, E) by numerical 
derivatives of P and M with respect to H and E, respec- 
tively (see Eqs. (1.1) or (2.8c)). In the inset of Fig. 3 it 
is shown that both data sets coincide. 

Figure 4 illustrates the temperature-driven melting of 
the singularity of the magnetoelectric tensor as function 
of the magnetic field. The values E = 1 oi the electric 
field and T ~ 0.001 J of the temperature belong to the 
red curve of the third panel of Fig. 3. With increasing 
temperature the curve becomes more and more flat. Note 
that the temperature in Fig. 4 is only slightly increased. 
For T = 0.03 J J, for instance, there is no longer any 
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FIG. 3: (Color online) Magnetization M, polarization P, and 
magnetoelectric tensor a as functions of the magnetic field 
H for anisotropy A = 4.5, temperature T/J = 0.001 and 
different values E = 0.1,1,2,3,4. In the inset of the third 
panel the two data sets JdnP and dsM are compared, both 
obtained numerically using the five-point stencil algorithm of 
the first derivative. 
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FIG. 4: (Color online) Temperature dependence of the mag- 
netoelectric tensor a for E = 1 and A = 4.5 close to the 
critical point H « 1.95J. The square-root singularity for 
T — decays rapidly with increasing temperature. 



indication of singular behavior. The two curves JdnP 
and OeM would overlay exactly such that only one of 
them is shown. 

In Figure 5 the temperature dependence of the elec- 
tric polarization and the magnetoelectric tensor are dis- 
played. The values of the electric and magnetic field 
arc chosen in a way such that the different ground-state 
phases are represented (see Fig. 1), just below and above 
the critical lines. All curves show non-monotonic be- 
havior and go to zero in the infinite-temperature limit. 
In the zcro-tcmpcrature limit the magnetoelectric tensor 
vanishes in the Neel-ordered and in the fully polarized 
ferromagnetic phase, whereas it remains constant in the 
Tomonaga-Luttingcr liquid. 

Figure 6 shows the dependence of the magnetization 
M and the polarization P on the electric field E for fixed 
values of the magnetic field H, again for very low temper- 
atures. For higher temperatures the cusps smooth out. 
Phase transitions only occur in the regimes (see Fig. 1 
with A = 4.5) 



H/J < Hi{E = 0)/J =2.637 . 
H/J > Hu{E = 0)/J = A + 1 



5.5. 



(3.9a) 
(3.9b) 



Here, the singularities are described by three different 
exponents. In most cases one can observe square-root 
behavior. This is always the case for the magnetiza- 
tion. However, for the special value of magnetic field 
H/J = 1 -|- A, which corresponds to the transition from 
the Tomonaga-Luttingcr liquid phase to the fully polar- 
ized ferromagnetic phase at zero electric field, the electric 
polarization grows quadratically in E. 

This can be explained as follows. The critical behavior 
of the correlation function (S'f S"! ) for fixed E as function 
of the magnetic field H is {SfS^) - ^/H~H for H < 
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FIG. 5: (Color online) Temperature dependence of the electric 
polarization P and the magnetoelectric tensor a close to the 
critical points H/J ^ 1.95 and H/J = 4.5 + \/2 (see Fig. 1). 



A + 1 and 



Hu- The gap between the fixed value H 
H^{E) = A + VT 

Due to the prefactor 



is quadratic in E. Therefore, the 
critical behavior is Hu 



H E 

in Eq. (3.5b) one eventually obtains 
2E 



P{E) 



E' 



ff=A + l 



H=A+i VTVW 

(3.10) 

In addition to that, the polarization as a function of 
the external electric field can also exhibit linear behav- 
ior with further transition to the square-root form, at 
the values of E corresponding to the transition from 
the Neel-ordered phase to the Tomonaga-Luttiger liquid 
phase. Thus, one can define the critical exponent 5e 
which characterizes the behavior of the polarization on 
its corresponding conjugate field (electric field), 



P ^\E~Er 



Se = 1/2,1,2. 



(3.11) 



The limiting value of P{E) ioi E oo can also be un- 
derstood. Due to Eqs. (2.4) and (2.6) the model reduces 
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FIG. 6: (Color online) Electric field dependence of the mag- 
netization M and the polarization P for different fixed values 
of the magetic field H at very low temperature T/J = 0.001. 
The values H/J = 2.637 and H/J = 5.5 correspond to the 
phase transitions from the Neel-ordered ground state to the 
Tomonaga-Luttinger liquid and from the Tomonaga-Luttinger 
liquid to the fully polarized ferromagnetic state, respectively 
(see Fig. 1). 



in this liniit to a free-fermion model with iJ — > 0. The 
rescaling J i-> J yields 



N 



(3.12) 



Then, the free energy per lattice site follows from 
Eq. (3.4a) and is given by 



I{T,H,E) 



T 

1" Jo 



los 



2 cosh 



/ Jcosfc\ 
V 2T j 



dfc. 



(3.13) 

Hence, the partition function Z = tr (e ^1^) ~ c 
with / = F/N and the nearest-neighbor correlation func- 
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tions can be calculated in the low-temperature limit, 
T 91og(Z) _ 1 df 



2N dJ 2 dJ 

1 r , /jcosfc\ , 

— / tann — — - — cos k ak 



4tt 



2T 



'^-^ r \cosk\dk^-^. (3.14) 



2tt 



Therefore, P{E) 1/tt ior E ^ oo and T ^ 0. 



IV. DISCUSSIONS AND CONCLUSION 

The system considered in the present paper, being an 
integrable one, allowed us to obtain an exact description 
of the MEE by means of QTM and NLIE techniques. 
We considered a microscopic dielectric polarization to be 
proportional to the antisymmetric product of x- and y- 
components of the spins from two adjacent sites, i. e. to 
a spin current . This mechanism, which is known as a 
KNB-mcchanism,^ is realized in several one-dimensional 
magnetic materials. ^^^^^ In our model, however, we have 
to impose a restriction in order to preserve integrability. 
While the standard model to describe the MEE in one- 
dimensional magnetic materials with KNB-mechanism is 
the Ji- J2-chain, Eq. (1.4), we consider an ordinary XXZ- 
chain in a magnetic field with DM-terms, Eq. (2.1). The 
DM-tcrms describe a local dielectric polarization where 
the corresponding coefficient is the electric field magni- 
tude E in appropriate units. 

The main objects of our investigation are the magneti- 
zation M and the dielectric polarization P, which are 
the thermal averages of the operators J2n=i ^.nd 
Y.n=iiS'nSn+i ~ S'nS'^_^_^), respectively. As a main re- 
sult the exact finite-temperature plots of M(T, H, E) and 
P{T, H, E) have been obtained, especially in the low- 
temperature regime. Below we summarize the most in- 
teresting features of these plots. 

As we mapped the system with DM-interaction onto 
the ordinary XXZ-chain, where E enters the effective 
coupling constant J{E) and effective anisotropy A(i?), 
the magnetization curves M{H) at fixed values of T and 
E (see Fig. 3, upper panel) do not show any great dif- 
ference from those of the conventional XXZ-chain. '^^ The 
polarization as a function of magnetic field at constant 
values of temperature and electric field exhibits non- 
monotonic behavior (see Fig. 3, middle panel). 

Since we chose a large anisotropy A ~ 4.5, one can 
distinguish three possible parts of the polarization curve 
corresponding to the Neel-ordered phase, Tomonaga- 
Luttinger massless regime, and fully polarized spin con- 
figuration of the effective XXZ-chain, respectively. All 
these three regimes are easily recognizable in the middle 
panel of Fig. 3. At low values of H, when the system is 
in the Neel-ordered ground state, the polarization P{H) 
is constant and thus starts with a plateau. One can also 



observe the effect of the electric field on the polariza- 
tion P{H). With increasing value of E the plateau be- 
comes shorter and higher. In the limit E oo the effec- 
tive anisotropy parameter A ~ l/vT+~E^ becomes zero, 
which means M — Q and P = Pmax = I/tt for iJ = 0. 
The second segment of the magnetic field dependent po- 
larization corresponds to the Tomonaga-Luttinger regime 
of the effective XXZ-chain. This part of the polarization 
curve is a dome which starts at the plateau and ends 
at the value P = corresponding to the fully polarized 
ground state of the spin system. The values of the crit- 
ical exponents corresponding to the two critical values 
of H (transitions from the Neel-ordered ground state to 
the Tomonaga-Luttinger regime and from the Tomonaga- 
Luttinger regime to the fully polarized spin state) are 
equal to 1/2. 

A feature of the model under consideration is the van- 
ishing polarization P{H) at E = 0. There only exists 
mutual influence of electric and magnetic properties of 
the model for E ^ 0. With respect to the electric prop- 
erties the system does not show ferroic order, but is a 
conventional 'dielectric' system which acquires polariza- 
tion in an external electric field. However, this dielectric 
polarization is strongly affected by the external magnetic 
field as well, which is the MEE one can observe in our 
system. 

Looking at the dependence of the magnetoelectric ten- 
sor a on the magnetic field H (see Fig. 3, lower panel) one 
can easily identify the zero-temperature phase transitions 
between Neel, Tomonaga-Luttinger and fully polarized 
ground states by the singular peaks. For sufficiently large 
H the polarization is zero due to the transition to the 
fully polarized ferromagnetic state. The magnetoelectric 
tensor has a negative peak at this transition point. In 
addition to that, the absolute magnitude of this peak is 
larger than that of the positive peak corresponding to 
the transition from the Neel-ordered to the Tomonaga- 
Luttinger liquid state. This indicates a more dramatic 
drop of polarization for upper critical magnetic fields 
compared to its rise for the lower critical fields. However, 
the type of singularity in both cases is the same with crit- 
ical exponent equal to 1/2, i.e. a ^ \H — Hc\ ■ 

The behavior of the polarization as function of elec- 
tric field is strongly affected by the magnetic field (see 
Fig. 6, lower panel). Even at vanishing magnetic field it 
exhibits non-trivial behavior during the whole process of 
polarization. The polarization increases linearly for small 
values of E followed by a smooth curve approaching its 
saturated value I/tt in the limit E -> oo. Switching on 
a magnetic field causes essential changes in the course of 
the polarization curve. For low magnetic fields the polar- 
ization starts linearly as in the H = case, but changes 
to a square-root behavior at the value of E correspond- 
ing to the transition from Neel-ordered ground state to 
Tomonaga-Luttinger liquid phase. At this point the po- 
larization curve has a cusp. For intermediate values of 
the magnetic field the system starts with a Tomonaga- 
Luttinger liquid phase, and the increase of E does not 
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drive the system to another ground state. One can ob- 
serve monotonic increase of the polarization with linear 
behavior at small E. At the special value of H corre- 
sponding to the critical value of the transition between 
Tomonaga-Luttinger liquid phase and fully polarized spin 
configuration a.t E = 0, i.e. H ^ + 1)^: one gets 
quadratic behavior for small E. Finally, for sufficiently 
large magnetic fields, when the system starts at the fully 
polarized magnetic configuration, the electric polariza- 
tion vanishes. This changes dramatically to a square- 
root increase at the point of quantum phase transition 
to Tomonaga-Luttinger liquid state. To summarize one 
can distinguish four different behaviors of the polariza- 
tion P{E): linear followed by a cusp for < H < Hi, 
linear with a smooth increase ior H ~ Q or Hi < H < Hu, 
quadratic for H — H^, and vanishing followed by a cusp 
for H > Hi. 

Finally, the temperature dependencies of the polariza- 
tion P and of the magnetoelectric tensor a illustrate an- 
other interesting feature (see Fig. 5). In the vicinity of 
quantum critical points the temperature dependence of 
the polarization is non-monotonic. The most intriguing 
feature shows the curve for H = 6J and E ~ 1, where at 
low temperatures polarization is absent (P = 0). Then, 
at intermediate temperatures the polarization curve rises 
up to values of about 0.065. Eventually, it drops to zero 
for high temperatures. Thus, starting from zero temper- 
ature and increasing it continuously, the system exhibits 
a broad polarization peak over almost five orders of mag- 
nitude. Here, the values of magnetic and electric fields 
correspond to the fully polarized ferromagnetic phase 
just above the critical line. The temperature dependence 
of the magnetoelectric tensor, which is presented in the 
lower panel of Fig. 5, is also non-monotonic and indicates 
the different regimes of the system response. 

Although the system considered in the present paper 
exhibits MEE and also demonstrates a number of inter- 
esting features like many regimes of polarization and a 
thermally activated peak, it does not exhibit ferroelectric 
polarization, which means that P is always zero, unless 
E ^ 0. In order to obtain more interesting and richer fea- 
tures like double ferroic order or induction of polarization 
only by the magnetization and/or by the magnetic field 
within the KNB-mechanism, one should consider more 
sophisticated models. One way which seems straightfor- 
ward to us is to include microscopic interactions between 
local order parameters, e. g. local magnetization and local 
polarization. One example is 

N 

Hint ^ ~ ■ (4-1) 

n=l 

The local polarization is composed according to KNB- 
mechanism and includes next-to-ncarcst-ncighbor spin 



interactions. Three-site interactions recently received 
considerable attention in a bit diverse context. It 
was shown by Suzuki in Ref. [55] that there is a series 
of spin chains with multiple-site interactions of special 
kind, which can be mapped onto free spinless fermions 
via the Jordan- Wigner transformation. Unfortunately, 
the S = 1/2 XXZ-chain with three-site spin- interaction 
terms given by Eq. (4.1) is no longer integrable. One has 
to restrict oneself to the XY-chain with corresponding 
terms^'^'^^^^^ or has to implement numerical simulations. 
Another obvious non-trivial choice of local interaction be- 
tween magnetization and polarization is X]n=i Pn'^^+i- 
However, in this case even the corresponding XX-model 
is not exactly solvable. 

Another possibility to construct the model exhibiting 
the double ferroic order and MEE within the integrable 
systems is to consider rather sophisticated spin chains 
with special three-spin interaction discussed in Ref. [60] . 
All these ideas conveyed above could be the guide line 
to further developments in the MEE and multiferroics of 
integrable systems. 
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Appendix A: Exact determination of correlation 
functions 

The functions (p, uj, and uj' that determine all static 
correlation functions of the XXZ-chain are defined in 
terms of solutions of linear and non-linear integral equa- 
tions. They were termed the physical part of the 
problem, since the physical parameters like tempera- 
ture or magnetic field enter solely through these func- 
tions. We provide their definitions only for the case 
A > 1. The definitions for the case — 1 < A < 1 can 
be found in Ref. [54]. 

First of all let us define a basic pair of auxiliary func- 
tions as the solution of the non-linear integral equations^'^ 
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logb(a;) ^ - 
log b{x) ^ — 



2Jsh{r]) 
T 

2Jsh(77) 
T 



H 



d{x) - + * log (1 + (x) -n * log (1 + b) (x) , 
d{x) + — + K * log (1 + b) (x) - n+ * log (1 + b) (x) , 



(Ala) 
(Alb) 



where we denote with * the convolution (/ * g){x) = and are given by 
i /^^/2 fi^ ~ y)9{y)^y- Here we introduced the param- 



etcr rj defined by A = ch(77). Eqs. (Al) are valid for all 
77 > meaning that A > 1. The driving term d is given 
by its Fourier series, 



k— — oo 



ch{r]k) 



(A2) 



Note that the physical parameters temperature T, mag- 
netic field H, and coupling J enter only through the driv- 
ing terms of Eqs. (Al) into our formulae. The kernels n 



k— — c 



^ — 7i\k\-\-2i'mkx 

2ch{r]k) 



K^{x) = k{x ±17] ) , 



(A3a) 
(A3b) 



where 77 =77 — 6 with an arbitrary small number e > 0. 

Except for the auxiliary functions b and b wc need two 
more pairs of functions gjf^^ and g'Jj^^ in order to define 
If, uj, and w'. Both pairs satisfy linear integral equations 
involving b and b, 



g+{x) = ~d{x -1^) + K^ \ {x) - K- * -J^^{x) , (A4a) 



T-(x) — n * T 

9^ t^^ .+ „ 9t 



Qu (a;) = ~d{x -iy) + K,* ' 1_i (x) - k+ * " (x) , 

1 + b 1 + ^' 

g'ti^) = ~c+{x -u) + K* ^ I {x) -K- * ^ 1_^ ix) + I * -I _i (x) - r (x) , 

1 + b 1 + 1 + b i + t, 

g'-Ax) = -c_ix-,y) + K* -^^{x) - K+ * -j^{x) + / * -^^{x) l+-^{x) . 

1 + b 1 + 1 + b 1 + 

The functions / and c± are again given by their Fourier series, 



(A4b) 
(A4c) 
(A4d) 



4ch^(r7fc) 2ch2(r7A:) ' 



where we set sign(O) = for the sign function. 

The functions ip{n), uj(fj.i, (12), and uj' {]j,i , (12) that determine the explicit form of the correlation functions of the 
XXZ-chain can be written as integrals involving b, b, gl^\ and s'jj^^- The function 

,H = +r (.#^-r#^)d., (A6a) 



27ry_,/2 Vl + fa(x)-i l + b(a;)- 

determines the magnetization Al (T, H) = —^ip{0) which is the only independent one-point function of the XXZ-chain. 
The function 



w(j^i, 1^2) = -4k(i/2 - 1^1) + Krf{iy2 -vi) -d-^ I -\ ^"L, I {V2) 

Vl+^ 1+b / 



(A6b) 
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also determines the energy per lattice site of the XXZ-chain, (sj_isj + Sj_iS^ + Asj_js|) = sh(77)aj(0, 0)/4. The 
function uj'{fii,fi2) is defined as 

ll>' {1^1,1^2) 



^1+ 



y + ^ 1 + b ) 

- C- * — — — Y(i^2) - C+ * —{V2) , 

1+0 1+b 



(A6c) 



r 



where we set v = —iv and Vj = ^ivj. The functions K^i 
and Lri are determined by 



as well as their linear counterparts for g^^'' and 
were solved iteratively in Fourier space utilizing the fast 
Fourier transformation algorithm. The derivatives of 
g^t'^ and g'Jj^^ with respect to /i, needed in the compu- 
tation of the derivative of lo' satisfy linear integral equa- 
tions as well, which were obtained as derivatives of the 
equations for g\^^ and g'^i^^ ■ Taking into account deriva- 
tives is particularly simple in Fourier space. 



K^{x) 
Ln[x) 



sh(277) 



2 sin(a:; + irj) sin(x — iri) 

i sin(2x) 
2 sin(.T + irj) sm(x — irj) 



(A7a) 
(A7b) 



For the calculation of the magnetization and the polar- 
ization, the non-linear integral equations for fa and fa 
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